Motzkin

Adapted from: (3.6) and (3.19) of [BPT12]

[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R. Semidefinite Optimization and Convex Algebraic Geometry. Society for Industrial and Applied Mathematics, 2012.

The first explicit example of nonnegative polynomial that is not a sum of squares was found by Motzkin in 1967. By the Arithmetic-geometric mean, $ \frac{x^4y^2 + x^2y^4 + 1}{3} \ge \sqrt[3]{x^4y^2 \cdot x^2y^4 \cdot 1} = x^2y^2 $ hence $ x^4y^2 + x^2y^4 + 1 - 3x^2y^2 \ge 0. $ The code belows construct the Motzkin polynomial using DynamicPolynomials.

using DynamicPolynomials
@polyvar x y
motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2

\[ 1 - 3x^{2}y^{2} + x^{2}y^{4} + x^{4}y^{2} \]

The Motzkin polynomial is nonnegative but is not a sum of squares as we can verify numerically as follows. We first need to pick an SDP solver, see here for a list of the available choices.

using SumOfSquares
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver)
@constraint(model, motzkin >= 0) # We constraint `motzkin` to be a sum of squares

optimize!(model)

We see that the problem is detected as infeasible...

termination_status(model)

INFEASIBLE::TerminationStatusCode = 2

... and that the dual solution is a certificate of the infeasibility of the problem.

dual_status(model)

INFEASIBILITY_CERTIFICATE::ResultStatusCode = 4

Even if the Motzkin polynomial is not a sum of squares, it can still be certified to be nonnegative using sums of squares. Indeed a polynomial is certified to be nonnegative if it is equal to a fraction of sums of squares. The Motzkin polynomial is equal to a fraction of sums of squares whose denominator is $x^2 + y^2$. This can be verified numerically as follows:

model = SOSModel(solver)
@constraint(model, (x^2 + y^2) * motzkin >= 0) # We constraint the `(x^2 + y^2) * motzkin` to be a sum of squares

optimize!(model)

Now the problem is declared feasible by the solver...

termination_status(model)

OPTIMAL::TerminationStatusCode = 1

... and the primal solution is a feasible point, hence it is a certificate of nonnegativity of the Motzkin polynomial.

primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

One may consider ourself lucky to have had the intuition that $x^2 + y^2$ would work as denominator. In fact, the search for the denominator can be carried out in parallel to the search of the numerator. In the example below, we search for a denominator with monomials of degrees from 0 to 2. If none is found, we can increase the maximum degree 2 to 4, 6, 8, ... This gives a hierarchy of programs to try in order to certify the nonnegativity of a polynomial by identifying it with a fraction of sum of squares polynomials. In the case of the Motzkin polynomial we now that degree 2 is enough since $x^2 + y^2$ works.

model = SOSModel(solver)
X = monomials([x, y], 0:2)

6-element DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}:
 1
 y
 x
 y²
 xy
 x²

We create a quadratic polynomial that is not necessarily a sum of squares since this is implied by the next constraint: deno >= 1.

@variable(model, deno, Poly(X))

\[ ({\_}_{1}) + ({\_}_{2})y + ({\_}_{3})x + ({\_}_{4})y^{2} + ({\_}_{5})xy + ({\_}_{6})x^{2} \]

We want the denominator polynomial to be strictly positive, this prevents the trivial solution deno = 0 for instance.

@constraint(model, deno >= 1)
@constraint(model, deno * motzkin >= 0)
optimize!(model)

termination_status(model)

primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We can check the denominator found by the program using JuMP.value

value(deno)

\[ 23.255419067239504 + 12.293357118019848y^{2} + 12.293357283456736x^{2} \]

Because a picture is worth a thousand words let's plot the beast. We can easily extend Plots by adding a recipe to plot bivariate polynomials.

using RecipesBase
@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)
    x, y, (x, y) -> p(variables(p) => [x, y])
end
import Plots
Plots.plot(
    range(-2, stop=2, length=100),
    range(-2, stop=2, length=100),
    motzkin,
    st = [:surface],
    seriescolor=:heat,
    colorbar=:none,
    clims = (-10, 80)
)

This page was generated using Literate.jl.